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ABSTRACT 

We describe a numerical code to solve the equations for ideal magnetohydrodynamics 
(MHD). It is based on an explicit finite difference scheme on an Eulerian grid, called the 
Total Variation Diminishing (TVD) scheme, which is a second-order-accurate extension 
of the Roe-type upwind scheme. We also describe a nonlinear Riemann solver for ideal 
MHD, which includes rarefactions as well as shocks and produces exact solutions for 
two-dimensional magnetic field structures as well as for the three-dimensional ones. The 
numerical code and the Riemann solver have been used to test each other. 

Extensive tests encompassing all the possible ideal MHD structures with planar sym- 
metries (i.e. one-dimensional flows) are presented. These include those for which the 
field structure is two-dimensional (i.e., those flows often called "1 + 1/2 dimensional") 
as well as those for which the magnetic field plane rotates (i.e.,, those flows often called 
"1 + 1/2 + 1/2 dimensional"). Results indicate that the code can resolve strong fast, slow, 
and magnetosonic shocks within 2-4 cells while more cells are required if shocks become 
weak. With proper stiffening, rotational discontinuities are resolved within 3-5 cells. Con- 
tact discontinuities are also resolved within 3-5 cells with stiffening and 6-8 cells without 
stiffening, while the stiffening on contact discontinuities in some cases generates numerical 
oscillations. Tangential discontinuities spread over more than 10 cells. 

Our tests confirm that slow compound structures with two-dimensional magnetic field 
are composed of intermediate shocks (so called "2-4" intermediate shocks) followed by slow 
rarefaction waves. Finally, tests demonstrate that in two-dimensional magnetohydrody- 
namics fast compound structures, which are composed of intermediate shocks (so called 
"1-3" intermediate shocks) preceded by fast rarefaction waves, are also possible. 

subject headings: hydromagnetics - magnetohydrodynamics:MHD - methodsmumerical - 
shock waves 



3 



1. INTRODUCTION 

Many astronomical objects as diverse as planets, stars, and galaxies all possess mag- 
netic fields which have important implications on their dynamics and evolution. However, 
in many objects including the Earth the ohmic dissipation time tp = L/r/ 2 (rj is the electric 
resistivity) is smaller than their ages. So their magnetic fields should be continuously gen- 
erated by some dynamo activity, otherwise the observed strength cannot be maintained. 
The dynamo activity involves patterns of magnetohydrodynamic (MHD) flows (for exam- 
ple, the interaction of differential rotation and convection as in some stars or accretion 
disks) needed to produce spatially coherent magnetic fields of large scale. However, such 
flows are usually highly nonlinear, so with several exceptions analytic approaches to un- 
derstand dynamo theory and the origin of the magnetic field have failed (for review and 
further references, see Parker 1979). 

With the appearance of fast supercomputers it has become possible to study the MHD 
flows numerically, but the development of numerical techniques to solve MHD equations has 
been slower due to the intrinsic complexity of the MHD flows. For instance, until now most 
numerical schemes to solve compressible MHD equations have been based on the methods 
using artificial viscosity (e.g.DeVore 1991; Lind, Payne & Meier 1991; Stone et al. 1992; 
Stone & Norman 1992) while schemes to solve compressible hydrodynamic equations have 
been based on methods using more sophisticated linear or nonlinear Riemann solvers 
(e.g. Roe 1981; Colella & Woodward 1984). It is rather recent that Brio & Wu (1988; 
henceforth BW) and Zachary & Colella (1992) developed schemes for MHD equations 
based on linear Riemann solvers and Dai & Woodward (1994; henceforth DW) one based 
on a nonlinear Riemann solver. These new schemes have been proved to handle MHD 
flows better, even though they are more expensive in CPU time and more difficult to code. 

In this paper, we describe a one-dimensional code to solve numerically the ideal MHD 
equations. It is based on an explicit second-order-accurate finite difference scheme on an 
Eulerian grid, called the Total Variation Diminishing (TVD) scheme, which was originally 
developed for numerical hydrodynamics by Harten (1983). It utilizes a Roe-type linear 
Riemann solver described in BW. Here, we do not attempt to describe the theoretical 
philosophy of the scheme which can be found in Harten (1983) and BW. Instead, we focus 
on the description of the procedure necessary to write the code and the test results for its 
performance. 

We also describe a nonlinear Riemann solver for ideal magnetohydrodynamics (MHDs) 
which has been developed to provide analytic solutions for shock tube tests of the code. 
It has been built by following a procedure similar to that in DW, but using initial guesses 
with the linear eigenvectors from BW and Zachary & Colella (1992). 

The paper is organized in the following way. In §2 we describe the step by step 
procedure of the code. In §3 we describe the nonlinear Riemann solver in details. We 
intend to make the detailed descriptions of the code and the nonlinear Riemann solver so 
others can recover them by following the descriptions. The tests using shock tubes are 
presented in §4 to demonstrate the ability of the code to capture discontinuities as well 
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as to follow smooth flows. The tests also serve to show the performance of the nonlinear 
Riemann solver. Finally, discussion is followed in §5 including the possible existence of fast 
compound structures as well as slow compound structures for two-dimensional magnetic 
field structures. 



2. TVD CODE FOR IDEAL MHDS 



2.1. Ideal MED Equations 



The subject of MHDs describes the dynamics of electrically conducting fluids in the 
presence of magnetic fields. The MHD equations represent coupling of the equations of fluid 
dynamics with the Maxwell's equations of electrodynamics. By neglecting the displacement 
current, the separation between ions and electrons, and the effects of electrical resistivity, 
viscosity, and thermal conduction, we get the following ideal MHD equations: 



| + v. W = o, 

+ v - Vv + -Vp (V x B) x B = 0, 

dt p p 

dp 

— + v • Vp + 7pV ■ v — 0, 
ot 

dB 

V x (v x B) 0. 



dt 



(2.1) 
(2.2) 

(2.3) 
(2.4) 



with an additional constraint V ■ B = for the absence of magnetic monopole (for details, 
see Shu 1992). Here, we have chosen units so that factor of An does not appear in the 
equations. 



For plane-symmetric, or one-dimensional flows exhibiting variation along the x direc- 
tion. The equations in (2.1) to (2.4) can be written in conservative form as 



dq dF 

dt dx 



= 0, 



(2.5) 



q = 



( p \ 

pv x 

pVy 

pv z 

By 

B z 
\ E J 



(2.6) 
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pv 2 +p* 



\ 



B 



X 



BxBy 



pv x vy 

F = pv x v z - B X B Z 

ByV X - B x Vy 

Bz^x — B x vz 

\(E + P*)V X - B X (B X V X + ByVy + B z v z ) ) 

where the total pressure and the total energy are given by 

P*=P + l(B 2 x + B 2 y + B 2 z 



E = \p(vl + v 2 y + vl) + ^ + \ (B 2 + B 2 + B 2 ) . 



(2.7) 



(2.8) 



(2.9) 



With the state vector, q, and the flux function, F{q), the Jacobian matrix, A(q) = dF /dq, 
is formed. The above system of equations is called hyperbolic, since all the eigenvalues of 
the Jacobian matrix are real and the corresponding set of the right eigenvectors is complete 
(Jeffrey & Taniuti 1964). However, the eigenvalues may coincide in some limiting cases 
(BW). 



2.2. Eigenvalues and Eigenvectors for Plane- Symmetric MHD Equations 

The first step to build a code based on the Harten's TVD scheme (Harten 1983) for 
the hyperbolic system of equations in (2.5) is to get the eigenvalues and the right and left 
eigenvectors of the Jacobian A(q). The seven eigenvalues in nondecreasing order are 



ai = v x - Cf, 


(2.10) 


a 2 = v x - c a , 


(2.11) 


as = v x- c s , 


(2.12) 


a 4 = v x , 


(2.13) 


a^ = v x + c s , 


(2.14) 


a& = v x + c a , 


(2.15) 


a 7 = v x + Cf, 


(2.16) 



where ct, c a , c s are the fast, Alfven, and slow characteristic speeds. The above represent 
the seven speeds with which information is propagated locally by three MHD wave families 
and an entropy mode. The three characteristic speeds are expressed as 



(2.17) 
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Cf = 



1 I a 2 + B x + B y + B z + 



1 j 2 £?2 + £?2 + £?2 
2<" + 



1 



where a is the sound speed given by 



a z + 





The corresponding right eigenvectors are (see, e.g., Jeffrey & Taniuti 1964) 

/ 1 \ 

V X ± Cf 



R 



V x ±Cf 



VyT 



B X ByC J 



,2 

7" 



B x B 



v z T 



B 



yC f 



2 
/ 



f 



,2 

7" 



Vx+Vy+Vz , ,± 




=FS 2 sign(S a; 
±Sy sign (S^ 
5z 

VP 
By 

VP 

K =F (-B^lty - ByVz) s\gn(B x ) J 

( 1 \ 

B X ByC S 

Vy T P(ci-cl) 



Rv x ±c s 



B y c 



B x B z c s 



P(cj-c 2 a ) 
Bzcj 



v%+vl+v% 



+ hf J 
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R Vx — 



1 

v x 

Vy 

V Z 




,2 , „,2. 

'y 



(2.24) 



where 



± cj B x cj (B y vy + B z vz) 7 _2 



2 / 2 2\ 
T (c/-a j, 



7-1 



?(c|- C 2) 7 



(2.25) 
(2.26) 



Near the point where either = or B, 



0, the above set of the right 



eigenvectors is not well defined with columns becoming singular. By renormalizing the 
eigenvectors, the singularities can be removed (BW). The renormalized eigenvectors are 



Rv x ±cj — 



a j (v x ±Cf 
oifVy =F a s (3yc a sign(S 
ajv z =F a s /3 z ca sign(S. 
a s f3yCj 
VP 

VP 



\ 



x) 



X) 



a 



vl+vl+v 2 z 



+ 9 



± 



f J 



(2.27) 



R 



V x ±C a - 








TPz sign(B x 
±f3 y sign(S x 

§z_ 
VP 

VP 

V =F (PzVy - PyVz) sign(S x ) J 



(2.28) 
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Rv x ±Cs — 



Ot S 

ot s (v x ± c s ) 
a s v y ±aj(3 y a sign(B x ) 

a s v z ±aj(3 z a sign(5 x ) 



c f VP 
af(3 z a 2 



ifVP 



V as f +9s J 



( 



R Vx — 



1 

Vy 

v z 



V 2 X +Vy+V 2 Z 

2 



where 



i I f ( \ 7 — ^ ( 2 2\ 
9f = — ± ajCfV x =F a s c a sign(S a; ) (/^i^ + (3 z v z j H — \ c f ~ a J ' 
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7-2 

7 



# s = ^— j- ± oi s c s v x ±atfa sign(S a; ) + (3 z v z j + ^—^ a s [c s - a 



7 

Here a's and /?'s are given by 



7-2 



a f = 

a s = 
A/ = 



cj — a 2 



By 
B z 



At the points where B y = B z = 0, /?'s are defined as the limiting values, i.e., 



By = B z = 0. 



(2.29) 



(2.30) 



(2.31) 
(2.32) 

(2.33) 

(2.34) 
(2.35) 
(2.36) 



(2.37) 



Similarly, at the point where B y = B z = and B x /p = a 2 , a's are defined as 



B z 

ctj = a s = 1 if By = B z = and — = a z 



(2.38) 



Then the left eigenvectors, which are orthonormal to the right eigenvectors, h\ 
_ / (1) (2) ,(3) ,(4) (5) (6) ,(7) 

^v x ±c f — y Vx ±Cfl l V x ±Cfl l V x ±Cfl l V x ±Cfl l V x ±Cfl l V x ±Cf l V x ±Cf 



l(l) 



1 CKf 9 o 1 

*-cr?r =F — 

M + 2 



-^-av x sign(S x ) - — c s \j3 y v y + f3 z v z ) 



,(2) 



1 ctf 9 la/- 
■- — -a v x ± - — J -a sign (B x , 



(3) _ _l«f 2 1^£, 



02 2 



,(6) / 2 , 2 - 7 2 \ /— 



1 CK 



/„2 



_ ,(2) ,(3) ,(4) ,(5) ,(6) ,(7) \ 

^V x ±C a yv x ±C a ' b V x ±C a 1 h V x ±C a 1 l V x ±C a 1 h V x ±C a 1 b V x ±C a > b V x ±C a ) 



1^ n 

l v x ±c a - 



/( 3 ) ^ z ■ in \ 

l v x ±c a = Tysign^), 

l { Sc a = ±1^(5,), 

(5) _ gyp 



i(6) 

l v x ±c a 



I 



(7) 

v x ±c a 



0, 



■ -Rm — 

(2.39) 
(2.40) 
(2.41) 
(2.42) 
(2.43) 
(2.44) 
(2.45) 
(2.46) 
(2.47) 

(2.48) 
(2.49) 
(2.50) 

(2.51) 

(2.52) 

(2.53) 
(2.54) 
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where 



_ / (1) (2) ,(3) ,(4) (5) ,(6) (7) 

^V X ±C S I l V x ±C s 1 l V x ±C s > l V x ±C s > 'l^iCg 5 6 t> x ±C s 5 6 t> x ±C s 5 6 t> x = 



,(1) 1 «s 2 2 _,_ 1 



^2 



±c. 



— c a i> x sign(Sa;) + -J-Cf(P y Vy + p z v z ) 



I 



(2) l«s 2 • • (-n . 

vJ± Cs = -^Y c f Vx ± ^T Ca Slgn ^' 



1 a s 

r 

1 a f 



,(3) 1 «s 2 ■ / i 



,(4) 



1 a s 



1 CK 



9 1 2 



02 2 



(5) 

i; x ±c s 


1 ttf 




(6) 

v x ±c s 


1 Ott 






,(?) _ 
l v x ±c s 


1 a s 2 



_ / (1) (2) ,(3) (4) (5) (6) (7) 
-^^x — I l v x , <-v x , lv x i <-v x , <-v x , h x , h x 



^ i a 2 a 2 + a 2 c 2 

lv x = 1 



e 



1 



(21 1/99 o o\ 
i (3) 1 / 2 2, 2 2\ 

^ = 77" («/ a + a s c /J^ 

r( 4 ) 1 / 2 2, 2 2\ 
* = ^ \ a .f a + a sCf)V Zl 

l vJ = ^a f a s pyCf (cj - c 2 ) y/p, 



4? = ^ a f a sPzCf (cj - C 2 S 
i (7) 1/22, 22 



/) 2 2/ 2, 2 -7 2^ , 22/2,2-7 2^ 

6>i = a /a + —0 J + a s c f {c s + —a ^ 

6> 2 = ajcja sign(Sa;) + a s c s c a sign(B x ), 
2 2,2,2 



(2.55) 
(2.56) 
(2.57) 
(2.58) 
(2.59) 
(2.60) 
(2.61) 
(2.62) 
(2.63) 

(2.64) 
(2.65) 
(2.66) 
(2.67) 
(2.68) 
(2.69) 
(2.70) 

(2.71) 

(2.72) 
(2.73) 
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Note that, with the above normalization, some columns are not continuous. In order 
to force them to be continuous, the following term 



sign(B T ) = 



1, 



if By > or By 



and B z > 



-1, if By < or By = and B z < 



(2.74) 



is multiplied R Vx ±c s and L Vx ± Ca if a 2 > c 2 and to R Vx ± Cf and L Vx ± Cf if a 2 < c 2 . 



It is interesting to see how the eigenvectors for the MHD equations reduce into those for 

I 1 1 \ 

the hydrodynamic equations in the limit B y — > and B z — > 0. If a > c^, R Vx ±cf ~ * Ry x ±a 
and R Vx — > -R^ D , that is, the characteristics associated with MHD fast waves become 
those associated with hydrodynamic sound waves. The other two eigenvectors for the 
three-dimensional hydrodynamic equations (Roe 1981) are obtained from the following 
combinations of those for slow and Alfven waves, 



V2 sign(S x ) 



R 



Vx+Cs 



— Rn„.— 



Vx-c s 



- (R Vx+Ca - Rvx-Ca) 



I o \ 



1 





\ v yJ 



(2.75) 



V2 sign(S x ) 



R 



Vx+Cs 



R 



Vx-Cs 



a 



+ (R Vx+Ca - R Vx - Ca ) 



Similarly, if a z < c^, R Vx ±c s —> R^± a and R Vx — > R^ x u , that is, the characteristics 
associated with MHD slow waves become those associated with hydrodynamic sound waves. 
The other two eigenvectors for the three-dimensional hydrodynamic equations are obtained 
from the following combinations of those for fast and Alfven waves, 



HD 



/ \ 


1 



\vz/ 



(2.76) 



V2 sign(S x ) 



R-vx+cj — R 



Vx-C 



f 



-(R 



Vx+Ca 



-R 



Vx-Ccu 



( o \ 

1 




\ v vJ 



(2.77) 



y/2 sign(S x ) 



R-Vx+Cj - Rvx-Cj 



+ (R 



Vx+Ca 



R 



Vx—Ca,/ 



( \ 



1 




\vz/ 



(2.78) 
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In the cases with purely two-dimensional magnetic fields and motions (i.e., B z = v z = 
0), the eigenvectors associated with Alfven waves becomes 



Rv x ±C a - 







±s\gn(B x ) 


l_ 

VP 
o 



\ 



(2.79) 



By combining them properly, they become (0, 0, 0, 1, 0, 0, 0) and (0, 0, 0, 0, 0, 1, 0) which 
are trivial. This indicates that the flows with two-dimensional magnetic fields do not 
produce structures associated with Alfven modes like rotational discontinuities (see §4 and 
5 for more discussions) 



2.3. A Second- Order Explicit TVD Scheme 



Here, we describe briefly the procedure to build the MHD-TVD code with the eigenval- 
ues and eigenvectors in the previous subsection. The purpose of this section is to provide 
a short but complete description of steps needed to build a code by the TVD scheme. 
For the details, e.g., why and how each step works, the choices for the values of internal 
parameters, etc., refer to the original work (Harten 1983). 

The convention for indices used in this subsection is the following. The superscript n 
represents the time step. The subscript i indicates the quantities defined in the cell centers 
while i + 1/2 identifies those defined on the cell boundaries. The subscript k represents the 
characteristic fields with the order that k = 1 is for the field associated with the eigenvalue 
v x — Cf, k = 2 for the field with v x — c a , k = 3 for the field with v x — c s , k = 4 for the field 
with v x , k = 5 for the field with v x + c s , k = 6 for the field with v x + c a , and k = 7 for the 
field with v x + Cf. 

In a code based on the TVD scheme, the physical quantities are defined in the cell 
centers while the fluxes are computed on the cell boundaries. Implementation of Roe's 
linearization technique would result in a particular form of the averaged physical quantities 
in the cell boundaries (Roe 1981). However, as pointed by BW, it is not possible to derive 
the analytic form of the averaged quantities in MHDs for general cases with the adiabatic 
index 7^2. Instead, we should modify the Roe's scheme by using the following simple 
averaging scheme, 

= P ~^p^, (2.80) 



v y,i+h 
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+ Vy,i+1 




2 


v z ,i 


+ v Zt i+i 




2 








2 


B z ,i 


+ B Z}i+ i 



Z,l+T 



2 



(2.82) 
(2.83) 
(2.84) 
(2.85) 
(2.86) 



= Pi +Pi+l 

2 2 

Then, other quantities like momentum, gas pressure, total energy, etc are calculated by 
combining the above quantities. Our tests for the cases with 7 = 2 indicated that the 
above simple averaging would do just as well when compared to the full implementation 
of Roe's linearization technique. 



The state vector q at the cell center is updated by calculating the modified flux / at 
the cell boundaries as follows: 



,n+l _ „n 



Qi - 



At" 
Ax 



fi+h = o + *W+l) 



'At 



(/i+± - fi-i), 



Ax ^ 

k=l 



9k,i+i 9k,i 



a 



0, 



for a 
for a 



= 



^ = sign(^. + i) max 0, min^|^. + i|, ^isign^+i) 



~ 2 



2 



l X, for|x|>2£_ 



£ = < 



r 0.1, for Jfe = 1 and 7 

0.2, for k = 2 and 6 

0.1, for /c = 3 and 5 

I. 0.0 for k = 4 



2.87) 

2.88) 

2.89) 
2.90) 

2.91) 

2.92) 
2.93) 

2.94) 
2.95) 
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Here, the time step At n is restricted by the usual Courant condition for the stability, 
At n = CcourAx/Maxflu" | + c" J with C COU r < 1. Typically we use C cour = 0.8. 

X,fc-|- 2 J ,A~r 2 

The rotational discontinuities represented by the k = 2 and 6 fields are steepened by 
replacing ^ , with g k l + d kji g kji : 



1 k,i 



k k — i I 



o, 



for d«A;,z+ll + l a A;,z-il) = 



(2.96) 



^ = sign^.+ i) max 0, min j sign(a fcji+ i)<T fc>i _ ia fe)i _i , tr^+ila^+i 



r fc,z+i 2 



1 



(2.97) 
(2.98) 



With the above steepening the rotational discontinuities are resolved within 3-5 cells, 
otherwise within 6-8 cells. The contact discontinuities represented by the k = 4 field 
could be also steepened by a similar scheme. However, our tests showed that the price for 
steepening the contact discontinuities is additional numerical oscillations. 



3. A NONLINEAR MHD RIEMANN SOLVER 



As a means to test quantitatively the MHD-TVD code and to understand more fully 
the properties of ideal MHD flows along one-dimension, we developed an accurate nonlinear 
MHD Riemann solver. It is similar in many respects to that described in DW. The 
most important improvement in our MHD Riemann solver is that it treats fast and slow 
rarefactions properly, instead of approximating them as "rarefaction shocks" as did the 
solver presented in DW. We present here only essential details of the "solver", referring 
readers to DW for other relevant background. 

As with hydrodynamic Riemann solvers, the construction of an MHD Riemann solver is 
based on the idea that two adjacent arbitrary states will evolve into a set of uniform states 
separated by left and right facing shocks and rarefactions. For the MHD problems, there 
are a total of eight states including the original pair. They are separated by six structures 
representing left and right propagating shocks or rarefactions of the three wave families 
and a structure representing the entropy mode (or the characteristic field associated with 
the eigenvalues v x in the discussion of the previous section) . The initial boundary moves 
with the structure of the entropy mode which becomes a contact discontinuity or, in 
degenerate cases, a tangential discontinuity (see §4 for more discussion). In fact the 
eigenstates developed in §2 provide approximate solutions to this problem. The particular 
difficulty in the MHD Riemann problem is that the equations are not strictly hyperbolic nor 
strictly convex (see the discussion in BW). In practice this means that the wave speeds 
of two families may sometimes coincide, and that compound wave structures involving 
both shocks and rarefactions may sometimes develop. For the most part we can ignore 
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this aspect, and consider those situations as special cases. Otherwise we can approach the 
Riemann problem in pretty much the same manner as for the hydrodynamical problem, 
except for the larger number of waves to consider. 

The solution of the problem is obtained by finding the required set of fast and slow 
shock jumps and rarefactions together with rotational discontinuities that self-consistently 
lead to a proper contact discontinuity or tangential discontinuity at the "center" of the 
structure. Our procedure consists of taking an initial guess for the resolved states of the 
(six interior) zones, and then iterating towards the proper jump conditions of the contact 
discontinuity or the tangential discontinuity. For this discussion, the left initial state is 
identified with zone 1 and the right initial state is identified with zone 8. Fast structures 
separate zones 1 and 2 as well as zones 7 and 8. Rotational discontinuities separate zones 
2 and 3 as well as zones 6 and 7. Zones 3 and 4 and zones 5 and 6 are each separated 
by slow structures. The contact discontinuity or the tangential discontinuity demarcates 
zones 4 and 5. 



Discontinuities are easier to handle than rarefactions, so we outline the methods to 
compute discontinuities first. The problem is conceptually simpler in Lagrangian mass 
coordinates, where the jump conditions across a discontinuity are (e.g., equation (2.7), 
DW): 



w[v] = -H, 


(3.1) 


W[v x ] = \p* - B*], 


(3.2) 


W[Vy] = —B X [By], 


(3.3) 


W[v z ] = -B X [B Z ], 


(3.4) 


W[VB y ] = -B x [v y ], 


(3.5) 


W[VB Z ] = -B x [v z ], 


(3.6) 


[V x p*] - B X [B X V X + ByVy + B z v z ], 


(3.7) 



where [Q] = — Q u is the difference between the downstream and upstream values of 
a quantity Q, W = —(pv x ) u = —(pv x )d is the Lagrangian speed of the discontinuity in 
the mass coordinates, and V = 1/p. The other quantities are as defined in §2. Thus, 
given the upstream state and an estimate of the Lagrangian speed of the discontinuity, the 
downstream state can be computed. 

The fast shock speed Wj and the slow shock speed W s in the Lagrangian mass 
coordinates are given in DW as 



W% = -^-^ {(C 2 S + C) + Si) ± ^(Cf + C) + Stf - 4(1 + S )(ClCj - S 2 )] , (3.8) 



where Cj = pcf and C s = pc s . The upper (lower) sign refers to fast (slow) waves. The 
coefficients So, Si, £2 can be written in terms of the jump in tangential magnetic field 
across the shock, [Bj_], as 

So = -^(7-1)^, (3-9) 
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Si = \ {-(7 - 2)Ci^l + 2C G 2 - ( 7 - 4)Ci - 2 7 Q 2 } (3.10) 
, 2 = 1 iCMJf + (7 + 2)C,Cg[BJ + ^ + + + ^ _ ^ [HJ 



(3.11) 

where C = pa is the Lagrangian sound speed, C a = pc a is the Lagrangian Alfven speed, 
C_L = y/pBj_, and = B^ + B\. All quantities except [BjJ in equations (3.9) to 
(3.11) are referred to the state upstream of the shock. The expressions used in equations 
(3.9) to (3.11) are equivalent to those given by DW, but are somewhat simpler and in 
our experience have provided more robust behavior in the Riemann solver, particularly in 
switch-on or switch-off shocks. An expression analogous to equation (3.8) appropriate to 
the special case B x = (magnetosonic shocks) is also given by DW. As is well known, fast 
and slow shocks do not alter the plane of the magnetic field. 

Rotational discontinuities, which are not compressive, can be handled by setting the 
jump [v x ] = in equations (3.1) and (3.2). That leads necessarily to W a = ±C a . The 
jump conditions required at the contact discontinuity can be found by setting W = in 
equations (3.1) to (3.7). With these results, it is clear that if the jumps of [Bj_] across 
shocks and the rotations of [Bj_] across rotational discontinuities are known then it should 
be possible to exactly determine the structure of any ideal MHD Riemann problem that 
involves only discontinuous interfaces. 

DW included rarefactions in their MHD Riemann solver by assuming they could also 
be treated as discontinuities (i.e., as "rarefaction shocks"). So long as the rarefactions 
are weak this is reasonably accurate, but not exact. In fact, it is also straightforward 
to include fast and slow rarefactions exactly, just as for hydrodynamics. By conserving 
all Riemann invariants through the rarefactions except that associated with the particular 
wave involved, one can derive a simple set of differential equations to be integrated through 
the rarefactions. The transitions computed in this way then replace the jumps given in 
equations (3.1) to (3.7). 

The relations appropriate to right facing (upper sign) and left facing (lower sign) fast 
rarefactions are (e.g., Jeffrey 1966) 

w 7 + 1,- C ± C 2 S 7 + l g|(gg - Cl) 1 + 1 PP > 

C °~ 2 ^ P C (C*-Cl)- 2 VP C 2 C±Co ~ 2 C ' (3 - 12) 

u' = t 1 Cl - Cl = ± C " - Cl - ± JL_9^L9k (3 13) 
x y/pCf(c$ - cl) ± ^pc ± c f ± 1 + iclc fP > [6A6) 

U!,J - "' : "=F^=^. (3-14) 



cos ip sin ip yfp C f 
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For right and left facing slow rarefactions one finds 



7 + 1 C±C j _ 7 + 1 r C % C< s ~ C a) _ 7 + 1 



x T vpc s (cj-ci) ± ^pc ± c s ± 1 + ic}c sP ' [6Ab) 
u 'y - u * _ T 1 Ca _ (3.17) 

cos ip sin t/> y/p C s 

In these expressions, tan^ = B z /B y , and primes represent derivatives with respect to Bj_. 
It is also useful to have the relationship 

(C) - C 2 a )(C 2 s - C 2 ) = -Clc\. (3.18) 

Our procedure for obtaining an accurate solution to a Riemann problem is similar to 
that employed by DW. We utilize the following steps: 



1) As an initial guess, we found it convenient and mostly reliable to use the eigenvectors 
defined in §2. In particular if the state vector, q, defined in equation (2.6) is represented 
in each of the intermediate regions by q(i) (i = 2, 3, ...7) then an estimate of q(i) is simply 

7 

q(x) = q(8)- £ a k R k , (3.19) 
k=i+l 

where a k is defined as in equation (2.90), and q(8) is the initial right state. Occasionally 
the states found by equation (3.19) can have unphysical properties (i.e., jumps in Bj_ that 
violate equations (3.1) to (3.7)), so some simple physical constraints need to be applied. 



2) The quantities Bj_(2), Bj_(4), Bj_(7) and i/j(3) from that solution are applied to deter- 
mine jumps across each of the six waves (e.g. [Bj_] 1^2)5 subject to the constraints that 
Bj_(b) = Bj_(4) and ip(4 — 6) = -i/>(3). This is accomplished by starting with the two fast 
waves, followed by the two rotations of Sj_ and finally the two slow waves. 



3) For the solution to be considered exact, we require that the resulting jump conditions 
at the contact discontinuity precisely satisfy those expected from equations (3.1) to (3.7). 
In particular, we test if all the jumps [v x ], [vy], [v z ] and [p*] = 0. 

4) If the contact discontinuity jump conditions are not satisfied, we vary Bj_(2), Sj_(4), 
Bj_(7) and ip(3) using a Newton- Raphson scheme based on a numerical approximation 
to the associated Jacobian matrix for an improved guess in the quantities Bj_(2), Bj_(4), 
B ± (7) and ip(S). 
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5) The procedure, beginning with step 2), is then repeated until convergence is obtained. 
The accuracy of the initial guess is the single most important aspect controlling the num- 
ber of iterations required. Once a reasonably approximate solution is found, usually only 
a couple of iterations lead to very good convergence. 

When B x = the same scheme is applied, but in somewhat simplified form, since there 
are no slow wave features or rotational discontinuities and the constraint ip(A) = ip(5) is 
removed. 

With these procedures, we are able in most cases to obtain solutions such that the 
contact discontinuity jump conditions are very well satisfied, even to near the limits of 
machine accuracy. In practice, we establish a convergence criterion of 1CP 6 relative to the 
larger of a zone averaged v x or fast wave speed for velocities or total pressure for [p*]. Only 
for Riemann problems involving switch-off or switch-on waves, we are significantly limited, 
since an exact switch-on or switch-off feature requires either the upstream or downstream 
flow speed to exactly equal the Alfven speed. In those cases, even very small errors in 
[Bj_] lead to significant errors in W, and so the remaining downstream state variables. We 
found that we could obtain precise solutions only if we permitted the smaller next to 
such a wave to have a value ~ 10 of the larger Bj_. These waves are for almost any 
practical purpose indistinguishable from exact switch-waves, however. 

We will discuss tests of the Riemann solver along with the tests of the MHD TVD 
code in the next section, since we test them against each other in many cases. Suffice 
it to say here that, to start with, we examined our Riemann solver against all those 
solutions presented in DW and against a number of other examples kindly given to us by 
Dr. Dai which were also generated with the Riemann solver described in DW. In the cases 
involving only discontinuities, our solutions agree exactly with the DW results. When 
there are rarefactions, we find some differences. But those are attributable to the fact 
that we treated rarefactions exactly, whereas they did not. Mostly those differences are 
relatively small, indicating that treating rarefactions as shocks would produce a reasonable 
result in an approximate MHD Riemann solver. 

4. NUMERICAL TESTS 

To test the code described in §2 as well as the Riemann solver described in §3, we chose 
MHD shock tube problems including those considered in BW and DW. In all the tests, 
we set the adiabatic index 7 = 5/3 and used a one-dimensional box with x = [0, 1]. The 
numerical calculations were done with a Courant constant 0.8 and without the stiffening 
of the contact discontinuity. The results of the numerical calculations with the code are 
plotted as dots and the analytic solutions of the Riemann solver are plotted as lines. The 
plotted quantities are density, gas pressure, total (thermal, kinetic, and magnetic) energy, 
x- velocity (parallel to the direction of structure propagation), y- velocity, z- velocity, y- 
magnetic field, ^-magnetic field, and the orientation angle of tangential magnetic field 
(ip = tan _1 (i? z / By)) in the plane perpendicular to the propagation vector. Numerical 
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values of the analytic solutions in the regions between the structures (e.g. between the left 
moving fast shock and the left moving rotational discontinuity, etc) are listed in the tables 
with the same labels as the figures. To simplify the discussion of this section we refer to 
solutions as two-dimensional when the magnetic field remains in one plane through the 
entire structure, or three-dimensional when the field cannot be so described. 

The first set of tests (also found in DW) has been done with two-dimensional field and 
velocity structure in the x — y plane but without change in the direction of the tangential 
magnetic field (B y in these tests). Fig. la shows the solution of the MHD shock tube test 
with the left state (p, v x , v y , v z , B y , B z , E) = (1, 10, 0, 0, 5/^/4%, 0, 20) and the right 
state (1, -10, 0, 0, 5/V4n, 0, 1) with B x = 5/V4n at time t = 0.08. The plot shows a 
pair of fast shocks, a left facing slow rarefaction, a right facing slow shock and a contact 
discontinuity. Fig.l b shows the solution of the MHD shock tube test with the left state 
(p, v x , v y , v z , By, B z , E) = (1, 0, 0, 0, 5/V4n, 0, 1) and the right state (0.1, 0, 0, 0, 
2/^/4tt, 0, 10) with B x = 3/^/4tt at time t = 0.03. The plot shows one fast shock and one 
fast rarefaction, one slow shock and one slow rarefaction, and a contact discontinuity. As 
expected from the two-dimensional nature of the field and velocity structure, there is no 
rotational discontinuity in either of these flows. In the numerical calculations, fast shocks 
that are strong with a large parallel velocity jump, [v x ] , are resolved within 2-4 cells, while 
slow shocks that are weak with a small jump require more cells to be resolved. Contact 
discontinuity spread typically spread over 6-8 or so cells. 

The second set of tests (also from DW) involves three-dimensional field and velocity 
structure where the magnetic field plane rotates. The solution of the MHD shock tube 
test with the left state (p, v x , v y , v z , B y , B z , E) = (1.08, 1.2, 0.01, 0.5, Z.S/^fAn, 2/VHr, 
0.95) and the right state (1, 0, 0, 0, 4/Vfe, 2/ v / 4tt, 1) with B x = 2/ v / 4tt at time t = 0.2 is 
plotted in Fig. 2a. Fast shocks, rotational discontinuities, and slow shocks propagate from 
each side of the contact discontinuity. The solution of the MHD shock tube test with the 
left state (p, v x , v y , v z , By, B z , E) = (1, 0, 0, 0, 6/\^An, 0, 1) and the right state (0.1, 0, 
2, 1, 1/V4tt, 0, 10) with B x = 3/^/4n at time t = 0.035 is plotted in Fig. 2b. A fast shock, 
a rotational discontinuity, and a slow shock propagate from the left side of the contact 
discontinuity, while a fast rarefaction, a rotational discontinuity, and a slow rarefaction 
propagate to the right. The rotation across the initial discontinuity of the magnetic field 
generates two rotational discontinuities. As in the previous case, strong fast shocks are 
resolved within 2-4 cells, while weak slow shocks take more cells. With proper stiffening, 
the rotational discontinuities spread over only 3-5 cells while the contact discontinuity 
spreads over more cells. 

In the third set, handling of magnetosonic structures with vanishing tangential flow 
velocity and parallel magnetic field is tested. A test from DW for magnetosonic shocks is 
set up with the left state (p, v x , v y , v z , B y , B z , E) = (0.1, 50, 0, 0, -1/^/4%, -2/\^4tv, 0.4) 
and the right state (0.1, 0, 0, 0, 1/V4tv, 2/ v / 4vr, 0.2) and with B x = 0. The solution has 
a pair of magnetosonic shocks propagating from a tangential discontinuity and is plotted 
at time t = 0.01 in Fig. 3a. The test for magnetosonic rarefactions are set up with the left 
state (p, v x , Vy, v z , B y , B z , E) = (1, —1, 0, 0, 1, 0, 1) and the right state (1, 1, 0, 0, 1, 0, 
1) and with B x = 0. The solution has only two identical magnetosonic rarefactions and is 
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plotted at time t = 0.1 in Fig. 3b. Magnetosonic shocks are fast shocks with zero parallel 
field (B x = 0) and they are resolved within 2-4 cells like fast shocks. The tangential 
discontinuity, which is a degenerate combination of one contact discontinuity, two slow 
structures, and two rotational discontinuities, spreads typically more than 10 cells. Even 
with the stiffening described in §2, we could not prevent the spreading of the tangential 
discontinuity in numerical calculations. But the Riemann solutions show that even slight 
differences from an exact tangential magnetic field, B x = 0, lead to the formation of the 
more complex and expanding set of features already mentioned. 

In the fourth set, we test how the MHD-TVD code and the Riemann solver deal with a 
special category of fast and slow structures; namely, the so-called "switch-on" and "switch- 
off" structures. The Tangential magnetic field turns on in the region behind switch-on fast 
shocks and switch-on slow rarefactions. On the other hand, it turns off in the region behind 
switch-off slow shocks and switch-off fast rarefactions. The test involving a switch-on fast 
shock has been set up with the left state (p, v x , v y , v z , B y , B z , E) = (1, 0, 0, 0, 1, 0, 
1) and the right state (0.2, 0, 0, 0, 0, 0, 0.1) and with B x = 1. The switch-on fast shock 
propagates to the right. Other structures formed in this test include a fast rarefaction, a 
slow rarefaction, a contact discontinuity, and a slow shock. The test results are plotted 
at time t = 0.15 in Fig. 4a, showing that the code and the Riemann solver handle the 
switch-on fast shock without any trouble. The test involving a switch-off fast rarefaction 
has been set up with the left state (p, v X: v y , v Z: B y , B z , E) = (0.4, —0.66991, 0.98263, 
0, 0.0025293, 0, 0.52467) and the right state (1, 0, 0, 0, 1, 0, 1) and with B x = 1.3. It has 
been designed to generate only a right-moving, switch-off fast rarefaction with a contact 
discontinuity with an accuracy (the ratio of residual strength to pre-rarefaction strength of 
tangential magnetic field) better than 0.3% (see Table 4b). The test results are plotted at 
time t = 0.15 in Fig. 4b. The numerical calculation shows small, yet noticeable, signatures 
of the right-moving slow structure and the left-moving hydrodynamic structure. The test 
involving a switch-off slow shock has been set up with the left state (p, v x , v y , v z , B y , 
B z , E) = (0.65, 0.667, -0.257, 0, 0.55, 0, 0.5) and the right state (1, 0.4, -0.94, 0, 0, 
0, 0.75) and with B x = 0.75. The solution of this test has, from left to right, a fast 
shock, a switch-off slow shock, a contact discontinuity, and a hydrodynamics shocks. In 
the region to the right of the switch-off slow shock hydrodynamic structures form with 
vanishing tangential magnetic field. The test results are plotted at time time t = 0.15 
in Fig. 4c, showing good agreement between the numerical calculation and the analytic 
solution. Our Riemann solver turns off the tangential magnetic field up to an accuracy 
of about 0.04% as indicated in Table 4c. The test involving a switch-on slow rarefaction 
has been set up with the left state (p, v x , v y , v z , B y , B z , E) = (1, 0, 0, 0, 0, 0, 1) and 
the right state (0.3, 0, 0, 1, 1, 0, 0.2) and with B x = 0.7. The structures formed in this 
test include a hydrodynamic rarefaction, a switch-on slow shock, a contact discontinuity, 
a slow shock, a rotational discontinuity, and a fast rarefaction. With the hydrodynamic 
nature in the region between the two left-moving rarefactions, no rotational discontinuity 
is produced on the left side of the contact discontinuity. Hence this test has a solution with 
only one rotational discontinuity, while two discontinuities are expected in most tests with 
three-dimensional field structure (e.g. the tests in Fig. 2). The test results are plotted at 
time t = 0.16 in Fig. 4d. 
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The final set of tests involves compound structures, whose existence was first discussed 
in BW. The initial setups includes two-dimensional field and velocity structure in the 
x — y plane and the sign of the tangential magnetic field (By) changes across the initial 
discontinuity. The test involving a slow compound structure has been set up with the left 
state (p, v x , Vy, v z , B y , B z , E) = (1, 0, 0, 0, 1, 0, 1) and the right state (0.125, 0, 0, 0, 
— 1, 0, 0.1) and with B x = 0.75. It is the same test as in BW, except we set the adiabatic 
index, 7 = 5/3, instead of 7 = 2. The results are plotted at time t = 0.1 in Fig. 5a. Here 
lines are the solution of the Riemann solver which has a rotational discontinuity and a slow 
shock instead of a slow compound structure, and numerical values of the solution are listed 
in Table 5a. Clearly the numerical calculations and the analytic solution produce different 
results around the compound structure, even though the agreement in other structures: 
two fast rarefactions, a slow shock, and a contact discontinuity, is acceptable for most 
purposes. Similarly, we set up a test involving a fast compound with the left state (p, 
v x , v y , v z , B y , B z , E) = (1, 0, 0, 0, 1, 0, 1) and the right state (0.4, 0, 0, 0, -1, 0, 0.4) 
and with B x = 1.3. The results are plotted at time t = 0.16 in Fig. 5b. Again lines are 
the solution of the Riemann solver with a rotational discontinuity and a fast rarefaction 
instead of a fast compound structure, and numerical values of the solution are listed in 
Table 5b. Further discussion on compound structures will follow in the next section. 

In the above tests with the MHD shock tube problems, the general and detailed 
agreement between the numerical calculations and the analytic solutions are satisfactory, 
and mostly excellent, except in the tests involving compound structures. Especially, the 
positions of shocks, rarefactions, and discontinuities agree very well, indicating the both 
generate the same solutions. Numerical results show that our code resolves strong shocks 
(fast, slow, or magnetosonic) typically within 2-4 cells, while more cells are required if 
shocks are weak with small jump in the parallel flow velocity. Rotational discontinuities 
are resolved within 3-5 cells with proper stiffening, and contact discontinuities are resolved 
within 6-8 cells without stiffening. While the stiffening of contact discontinuities make 
them look sharper in some cases, it generates spurious numerical oscillations in many cases. 
Tangential discontinuities spread over more than 10 cells. Further work could be done to 
develop proper stiffening schemes for contact discontinuities and tangential discontinuities 
to improve the ability of the MHD-TVD code to handle these discontinuities. 

5. DISCUSSION 

One test problem that has attracted considerable attention was one originally presented 
by BW involving a slow compound structure, which is related to structures known as 
"intermediate" shocks (e.g., Wu 1988; Kennel et al. 1990). This structure reverses the 
direction of and leads to a flow that passes from super-Alfvenic to sub-Alfvenic. 
Numerical MHD codes, which are always dissipative, seem generally to find the BW 
compound structure out of that shock-tube initial condition. For the compound wave 
it does not appear necessary that the magnetic field anywhere exist outside the upstream 
or downstream field plane. Thus, it would appear that this structure might be represented 
as a "1 + 1/2 dimensional flow problem". However, ideal MHD Riemann solvers that 
allow rotational discontinuities will generally find a solution to that shock-tube problem 
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that interprets the compound structure as a rotational discontinuity followed by a slow 
shock, making it a "1 + 1/2 + 1/2 dimensional" flow. 

It is not our intent to add to the controversy surrounding the reality of intermediate 
shocks. However, in order to test the performance as well as understand the properties of 
our code, we have done a high resolution calculation for the test of the slow compound 
structure in Fig. 5a with 8192 cells. The magnified region around the slow compound 
structure has been plotted in Fig. 6a. From the plot, the numerical values of several 
flow quantities on the both sides of the "shock" read: p = 0.647 (density), p g = 0.484 
(gas pressure), u = 0.963 (flow velocity in the shock frame), By = 0.536 (tangential 
magnetic field), Cj = 1.423 (fast speed), c a = 0.932 (Alfven speed), c s = 0.732 (slow 
speed), a = 1.117 (sound speed) in the preshock region, p = 0.827, p g = 0.758, u = 0.751, 
By = -0.0738, c f = 1.241, c a = 0.825, c s = 0.822, a = 1.236 in the postshock region, 
and shock velocity in lab frame is 0.284. In the preshock region the flow velocity is 
c a < u < cj (sub-fast but super- Alfvenic) with Alfvenic Mach number M a = 1.03, while 
in the postshock region u < c s (sub-slow) with M a = 0.910. Hence, the shock in the 
numerical calculation with our code should be considered as a "2-4" intermediate shock 
and the slow compound structure as a intermediate shock followed by a slow rarefaction. 

We have done also a high resolution calculation for the test of the fast compound 
structure in Fig. 5b with 8192 cells and the magnified region around the fast compound 
structure has been plotted in Fig. 6b. The numerical values in the plot of several flow 
quantities on the both sides of the "shock" read: p = 0.6820, p g = 0.5290, u = 1.617, 
By = 0.0555, cj- = 1.577, c a = 1.574, c s = 1.135, a = 1.137 in the preshock region, 
p = 0.7361, p g = 0.6031, u = 1.498, B y = -0.3442, cj = 1.622, c a = 1.515, c s = 1.092, 
a = 1.169 in the postshock region, and shock velocity in lab frame is 0.935. In the preshock 
region the flow velocity is u > ct (super-fast) with M a = 1.027, while in the postshock 
region c s < u < c a (sub-Alfvenic but super-slow) with M a = 0.989. Hence, the shock in 
the numerical calculation with our code should be considered as a "1-3" intermediate shock 
and the fast compound structure as a intermediate shock preceded by a fast rarefaction. 
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FIGURE CAPTIONS 

Fig. la. — Solution of the MHD shock tube test with the left state (p, v x , Vy, v z , By, B z , 
E) = (1, 10, 0, 0, 5/y/4n, 0, 20) and the right state (1, -10, 0, 0, 5/V&, 0, 1) with 
B x = 5/V47T and 7 = 5/3 at time t = 0.08 (test in DW table 7). Dots are the result 
of a numerical calculation with the MHD-TVD code described in §2 using 512 cells 
and a Courant constant of 0.8. Lines are the result with the nonlinear Riemann solver 
described in §3. Plots show from left to right (1) fast shock, (2) slow rarefaction, (3) 
contact discontinuity, (4) slow shock, and (5) fast shock. 

Fig. lb. — Solution of the MHD shock tube test with the left state (p, v x , Vy, v z , By, 
B z , E) = (1, 0, 0, 0, 5/v 7 ^, 0, 1) and the right state (0.1, 0, 0, 0, 2/ v / &, 0, 10) 
with B x = and 7 = 5/3 at time t = 0.03 (test in DW table 3a). Dots are the 

result of a numerical calculation with the MHD-TVD code described in §2 using 512 
cells and a Courant constant of 0.8. Lines are the result with the nonlinear Riemann 
solver described in §3. Plots show from left to right (1) fast shock, (2) slow shock, (3) 
contact discontinuity, (4) slow rarefaction, and (5) fast rarefaction. 

Fig. 2a. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, 
B z , E) = (1.08, 1.2, 0.01, 0.5, 2/v 7 ^, 0.95) and the right state (1, 0, 0, 

0, 4/^/4%, 1) with B x = 2/ v / 4tt and 7 = 5/3 at time t = 0.2 (test in DW 

table la). Dots are the result of a numerical calculation with the MHD-TVD code 
described in §2 using 512 cells and a Courant constant of 0.8. Lines are the result 
with the nonlinear Riemann solver described in §3. Plots show from left to right (1) 
fast shock, (2) rotational discontinuity, (3) slow shock, (4) contact discontinuity, (5) 
slow shock, (6) rotational discontinuity, and (7) fast shock. 

Fig. 2b. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, 
B z , E) = (1, 0, 0, 0, 6/v 7 ^, 0, 1) and the right state (0.1, 0, 2, 1, 1/V&, 0, 10) 
with B x = 3/^/4n and 7 = 5/3 at time t = 0.035 (test in DW table 5a). Dots 
are the result of a numerical calculation with the MHD-TVD code described in §2 
using 512 cells and a Courant constant of 0.8. Lines are the result with the nonlinear 
Riemann solver described in §3. Plots show from left to right (1) fast shock, (2) 
rotational discontinuity, (3) slow shock, (4) contact discontinuity, (5) slow rarefaction, 
(6) rotational discontinuity, and (7) fast rarefaction. 

Fig. 3a. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , B y , B z , 
E) = (0.1, 50, 0, 0, -l/v 7 ^, -2/v 7 ^, 0.4) and the right state (0.1, 0, 0, 0, 
2/ v / 4tt, 0.2) with B x = and 7 = 5/3 at time t = 0.01 (test in DW table 2a). Dots 
are the result of a numerical calculation with the MHD-TVD code described in §2 
using 512 cells and a Courant constant of 0.8. Lines are the result with the nonlinear 
Riemann solver described in §3. Plots show from left to right (1) magnetosonic shock, 
(2) tangential discontinuity, and (3) magnetosonic shock. 

Fig. 3b. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, 
B z , E) = (1, -1, 0, 0, 1, 0, 1) and the right state (1, 1, 0, 0, 1, 0, 1) with B x = 
and 7 = 5/3 at time t = 0.1. Dots are the result of a numerical calculation with the 
MHD-TVD code described in §2 using 512 cells and a Courant constant of 0.8. Lines 
are the result with the nonlinear Riemann solver described in §3. Plots show from left 
to right (1) magnetosonic rarefaction, and (2) magnetosonic rarefaction. 
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Fig. 4a. — Solution of the MHD shock tube test with the left state (p, v x , Vy, v z , By, B z , 
E) = (1, 0, 0, 0, 1, 0, 1) and the right state (0.2, 0, 0, 0, 0, 0, 0.1) with B x = 1 and 
7 = 5/3 at time t = 0.15. Dots are the result of a numerical calculation with the 
MHD-TVD code described in §2 using 512 cells and a Courant constant of 0.8. Lines 
are the result with the nonlinear Riemann solver described in §3. Plots show from left 
to right (1) fast rarefaction, (2) slow rarefaction, (3) contact discontinuity, (4) slow 
shock, and (5) switch-on fast shock. 

Fig. 4b. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, B z , 
E) = (0.4, -0.66991, 0.98263, 0, 0.0025293, 0, 0.52467) and the right state (1, 0, 0, 0, 
1, 0, 1) with B x = 1.3 and 7 = 5/3 at time t = 0.15. Dots are the result of a numerical 
calculation with the MHD-TVD code described in §2 using 512 cells and a Courant 
constant of 0.8. Lines are the result with the nonlinear Riemann solver described 
in §3. Plots show from left to right (1) contact discontinuity, and (2) switch-off fast 
rarefaction. 

Fig. 4c. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, B z , 
E) = (0.65, 0.667, -0.257, 0, 0.55, 0, 0.5) and the right state (1, 0.4, -0.94, 0, 0, 0, 
0.75) with B x = 0.75 and 7 = 5/3 at time t = 0.15. Dots are the result of a numerical 
calculation with the MHD-TVD code described in §2 using 512 cells and a Courant 
constant of 0.8. Lines are the result with the nonlinear Riemann solver described in 
§3. Plots show from left to right (1) fast shock, (2) switch-off slow shock, (3) contact 
discontinuity, and (4) hydrodynamic shock. 

Fig. 4d. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , B y , B z , 
E) = (1, 0, 0, 0, 0, 0, 1) and the right state (0.3, 0, 0, 1, 1, 0, 0.2) with B x = 0.7 
and 7 = 5/3 at time t = 0.16. Dots are the result of a numerical calculation with the 
MHD-TVD code described in §2 using 512 cells and a Courant constant of 0.8. Lines 
are the result with the nonlinear Riemann solver described in §3. Plots show from 
left to right (1) hydrodynamic rarefaction, (2) switch-on slow rarefaction, (3) contact 
discontinuity, (4) slow shock, (5) rotational discontinuity, and (6) fast rarefaction. 

Fig. 5a. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , By, B z , E) 
= (1, 0, 0, 0, 1, 0, 1) and the right state (0.125, 0, 0, 0, -1, 0, 0.1) with B x = 0.75 and 
7 = 5/3 at time t = 0.1 (test in BW). Dots are the result of a numerical calculation 
with the MHD-TVD code described in §2 using 512 cells and a Courant constant of 
0.8. Lines are the result with the nonlinear Riemann solver described in §3. Plots show 
from left to right (1) fast rarefaction, (2) slow compound, (3) contact discontinuity, 
(4) slow shock, and (5) fast rarefaction. 

Fig. 5b. — Solution of the MHD shock tube test with the left state (p, v x , v y , v z , B y , B z , 
E) = (1, 0, 0, 0, 1, 0, 1) and the right state (0.4, 0, 0, 0, -1, 0, 0.4) with B x = 1.3 
and 7 = 5/3 at time t = 0.16. Dots are the result of a numerical calculation with the 
MHD-TVD code described in §2 using 512 cells and a Courant constant of 0.8. Lines 
are the result with the nonlinear Riemann solver described in §3. Plots show from left 
to right (1) fast compound, (2) slow shock, (3) contact discontinuity, (4) slow shock, 
and (5) fast rarefaction. 

Fig. 6a. — Same slow compound structure in the MHD shock tube test as that in Fig. 5a. 
The calculation has been done with the MHD-TVD code described in §2 using 8192 
cells and only the region around the slow compound structure has been plotted. Plots 
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show that the slow compound structure from the numerical calculation is composed 
of a "2-4" intermediate shock followed by a slow rarefaction. 

6b. — Same fast compound structure in the MHD shock tube test as that in Fig. 5b. 
The calculation has been done with the MHD-TVD code described in §2 using 8192 
cells and only the region around the fast compound structure has been plotted. Plots 
show that the fast compound structure from the numerical calculation is composed of 
a "1-3" intermediate shock preceded by a fast rarefaction. 



